/*
 * CO2Solubility.cpp
 *
 *  Created on: 12 Sep 2011
 *      Author: Allan
 */

#include "CO2Solubility.h"

// C++ includes
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;

// GeoReact includes
#include "../Assembler.h"
#include "../Applications/Equilibrator.h"
#include "../Utilities/Range.h"

double CalculateCO2Solubility(double T, double P, double mNaCl);
void OutputCO2Solubility(double Tmin, double Tmax, double Tstep, const vector<double>& Ps, double mNaCl);

void DemoCO2Solubility()
{
	OutputCO2Solubility(30.0, 150.0, 1.0, {10.0, 50.0, 100.0}, 4.0);
}

double CalculateCO2Solubility(double T, double P, double mNaCl)
{
	System system = AssembleSystem({"H2O", "NaCl(aq)", "CO2(g)", "H2O(g)"});
	
	system.SetTemperature(T);
	system.SetPressure(P);
		
	Equilibrator equilibrator;
	
	equilibrator.ResetConditions();
	
	equilibrator.AddCondition("ChargeBalanced");
	equilibrator.AddCondition("FixedMass", "H2O", 55.508);
	equilibrator.AddCondition("FixedTotalMass", "Na[+]", mNaCl);
	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", mNaCl);
	equilibrator.AddCondition("FixedTotalMass", "HCO3[-]", 3.0);
	
	equilibrator.Equilibrate(system);
		
	return 3.0 - system.Get("Moles", "CO2(g)");
}

void OutputCO2Solubility(double Tmin, double Tmax, double Tstep, const vector<double>& Ps, double mNaCl)
{
	cout << setw(15) << left << "T";
	
	BOOST_FOREACH(double P, Ps)
		cout << "P = " << setw(11) << left << P;
	
	cout << endl;
	
	BOOST_FOREACH(double T, Range(Tmin, Tmax, Tstep))
	{
		cout << setw(15) << left << T;
		
		BOOST_FOREACH(double P, Ps)
			cout << setw(15) << left << CalculateCO2Solubility(T, P, mNaCl);
		
		cout << endl;
	}
}

